suppressPackageStartupMessages(library(dplyr))
library(readr)
library(tibble)
library(tidyr)
library(stringr)
library(forcats)
library(variancePartition)
library(BiocParallel)
library(edgeR)
library(fgsea)
library(readxl)
library(writexl)
library(ggplot2)
library(ggpubr)
#library(DOSE)
library(ggrepel)
library(furrr)
library(purrr)
library(future)
library(plotly)
set.seed(1969)
theme_set(theme_minimal())
register(MulticoreParam(workers=64))
plan(multisession, workers = 64)
dgea_label <- function(res) {
if (nrow(res) == 0) return(res)
res$diffexpressed <- "NO"
# if log2Foldchange > 1 and pvalue < 0.05, set as "UP"
res$diffexpressed[res$logFC > 1 & res$adj.P.Val < 0.05] <- "UP"
# if log2Foldchange < -1 and pvalue < 0.05, set as "DOWN"
res$diffexpressed[res$logFC < -1 & res$adj.P.Val < 0.05] <- "DOWN"
# Now write down the name of genes beside the points...
res$label <- NA_character_
res$label[res$diffexpressed != "NO"] <-
res$gene_name[res$diffexpressed != "NO"]
return(res)
}
plot_volcano <- function(tab, t) {
if (nrow(tab) > 0) {
ggplot(data=tab, aes(x=logFC, y=-log10(adj.P.Val), color = diffexpressed, label=label)) +
geom_point() +
scale_color_manual(values=c("#0006b1", "grey", "#bb0c00")) +
geom_vline(xintercept=c(-1, 1), col="grey", linetype = 'dashed') +
geom_hline(yintercept=-log10(0.05), col="grey", linetype = 'dashed') +
theme_minimal() +
geom_label_repel( max.overlaps = 30) +
ylab(expression("-log"[10]*"Adjusted p-value")) +
xlab(expression("log"[2]*"FC")) +
theme_pubr() +
theme(legend.position ="none") +
labs_pubr()+
labs(title = t) +
theme(plot.title = element_text(hjust = 0.4))
} else {
ggplot() +
annotate("text", x = 0.5, y = 0.5, size = 8,
label = sprintf("No significant genes for %s", t)) +
theme_void()
}
}
datapath <- "/sc/arion/projects/load/users/goateomics/garref01/garref01.2024-11-11.B2_B3_KD_apoepop_iMGL_bulkRNAseq"
dir.create("output", showWarnings = FALSE, recursive = TRUE)
dir.create("input", showWarnings = FALSE)
dir.create("shiny/www", showWarnings = FALSE)garref01.2024-11-11.B2_B3_KD_apoepop_iMGL_bulkRNAseqDream
Experimental design
B2 B3 KD apoe pop iMGL bulk RNASeq
96 samples total
- APOE: 4 lines
- 33 = ID1, ID3, ID4, ID6
- 44 = ID8, ID9, ID11, ID13
- siRNA: 3 genotypes
- SCR
- B2
- B3
- stimulus
- none
- THPG
- expt (two experiments/line)
- 1
- 2
Setup environment
Get protein coding genes from GENCODE GTF
Ensure that the gencode version of the GTF read in here matches that of the gtf used for STAR/Salmon, or things will not work properly.
read_gtf <- function(gtf_fname) {
rtracklayer::import(gtf_fname) |>
as_tibble() |>
select(type, gene_type, gene_id, gene_name, transcript_id, transcript_name)
}
gencode45_gtf <- read_gtf("data/gencode.v45.basic.annotation.gtf.gz")
get_protein_coding <- function(gtf_tab) {
gtf_tab |>
filter(type == "gene" & gene_type == "protein_coding") |>
select(gene_id, gene_name)
}
gencode45_genes <- get_protein_coding(gencode45_gtf)Read and prep data
Data are from: /sc/arion/projects/load/users/goateomics/garref01/garref01.2024-11-11.B2_B3_KD_apoepop_iMGL_bulkRNAseq
Count matrix (gene)
#import count metrix
cm_raw <- read_tsv("data/salmon.merged.gene_counts.tsv",
col_types = cols(.default = "d",
gene_id = "c",
gene_name = "c"))
# Ensure that the sample IDs in your sample sheet and the count matrix match
# modify this code to do so
cm <- cm_raw
head(cm)Sample sheet
# getting data from path in minerva
ss_raw <-
file.path(datapath, "deg_samplesheet.csv") |>
read_csv(col_types = cols(.default = "c", RIN = "d"))
ss_raw# Turn categorical variables into factors
proc_samples <- function(ss, bin_vars, grp, relevel) {
binary_vars <- unique(c(bin_vars, names(relevel)))
for (var in binary_vars) {
if (var %in% names(relevel)) {
ss[var] <- fct_relevel(ss[[var]], relevel[[var]])
} else {
ss[var] <- as_factor(ss[[var]])
}
}
for (group in grp) {
gname <- paste(group, collapse = "_")
ss <- ss |>
unite(!!gname, all_of(group), sep = ":", remove = FALSE) |>
mutate({{ gname }} := as_factor(!!sym(gname)))
}
out <- ss |>
mutate(sample = make.names(sample)) |>
column_to_rownames("sample")
ccols <- colnames(select_if(out, is.character))
if (length(ccols) > 0) {
warning(paste("The following columns are character, not factor: ", ccols))
}
return(out)
}
ss <- ss_raw |>
select(-APOE_siRNA, -APOE_siRNA_stim, -starts_with("fastq")) |>
proc_samples(bin_vars = c("APOE", "line", "siRNA", "stimulus", "expt"),
grp = list(c("APOE", "siRNA"), c("APOE", "siRNA", "stimulus"),
c("line", "expt")),
relevel = list(siRNA = "SCR", stimulus = "none", APOE = "33"))
ssgenenames <- cm |> select(gene_id, gene_name)
head(genenames)Exploratory data analysis
Make edgeR DGE list for exploratory analysis
Filter to expressed genes only
make_dgelist <- function(ss, cm, genes, grp, norm_fun = calcNormFactors) {
cm_mat <- cm |>
semi_join(genes, by = "gene_id") |>
column_to_rownames("gene_id") |>
select(rownames(ss)) |>
as.matrix()
#normalize DGE by reads (whatever that is calculated by) across samples
d_all <- DGEList(cm_mat)
d_all_norm <- norm_fun(d_all)
matdim <- list(samp = ncol(d_all))
matdim["all"] <- nrow(d_all)
#remove lowly expressed genes across samples, recalculate library size after removing genes
keep <- filterByExpr(d_all_norm, group = grp)
d <- d_all_norm[keep, , keep.lib.sizes = FALSE]
# should this be d <- d_all_norm[keep, , keep.lib.sizes = TRUE]? We need this if using CPM instead of just TMM
matdim["expr"] <- nrow(d)
return(list(d = d, ss = ss, matrix_dims = matdim))
}
d_all <- ss |>
make_dgelist(cm, gencode45_genes, "APOE_siRNA_stimulus")Print MDS Plots
extract_mds <- function(dmat, ss, n = 16) {
mds <- plotMDS(dmat, plot = FALSE)
eigenvec <- mds$eigen.vectors[seq_len(nrow(mds$eigen.vectors)), seq_len(n)]
eigenval <- mds$eigen.values[seq_len(n)]
md_tab <- {eigenvec %*% diag(sqrt(eigenval))} |>
as_tibble(.name_repair = \(n) paste0("MD", seq_along(n)))
md_tab_ss <- bind_cols(as_tibble(ss, rownames = "sample"), md_tab)
list(mds = md_tab_ss, var_explained = mds$var.explained[seq_len(n)])
}
mds <- extract_mds(d_all$d, d_all$ss)
plot_mds <- function(mds, grp) {
pct <- round(mds$var_explained * 100)[1:2]
p <- mds$mds |>
ggplot(aes(x = MD1, y = MD2, text = sample, color = !!sym(grp))) +
geom_point() +
labs(title = grp,
x = sprintf("Leading logFC dim 1 (%i%%)", pct[1]),
y = sprintf("Leading logFC dim 2 (%i%%)", pct[2]))
ggplotly(p, tooltip = "text")
}
plot_mds_34 <- function(mds, grp) {
pct <- round(mds$var_explained * 100)[3:4]
p <- mds$mds |>
ggplot(aes(x = MD3, y = MD4, text = sample, color = !!sym(grp))) +
geom_point() +
labs(title = grp,
x = sprintf("Leading logFC dim 3 (%i%%)", pct[3]),
y = sprintf("Leading logFC dim 4 (%i%%)", pct[4]))
ggplotly(p, tooltip = "text")
}plot_mds(mds, "APOE")plot_mds_34(mds, "APOE")plot_mds(mds, "siRNA")plot_mds_34(mds, "siRNA")plot_mds(mds, "line")plot_mds(mds, "stimulus")plot_mds(mds, "expt")Differential Gene Expression
Variance partition analysis
#given all these factors, where the model attributes variance to, everything you dont explain will end up in the residuals and will determine the SEM
f <- ~ (1 | APOE) + (1 | siRNA) + (1 | stimulus) + (1 | line) + (1 | expt)
vm_vp <- voomWithDreamWeights(d_all$d, f, d_all$ss, BPPARAM = bpparam())Memory usage to store result: >139.7 Mb
Dividing work into 128 chunks...
Total:9 s
vp <- fitExtractVarPartModel(vm_vp, f, d_all$ss, BPPARAM = bpparam())Dividing work into 128 chunks...
Total:9 s
saveRDS(vp, "input/vp.rds")
plotVarPart(sortCols(vp))#given all these factors, where the model attributes variance to, everything you dont explain will end up in the residuals and will determine the SEM
f <- ~ (1 | APOE) + (1 | siRNA) + (1 | stimulus) + (1 | line_expt)
vm_vp <- voomWithDreamWeights(d_all$d, f, d_all$ss, BPPARAM = bpparam())Memory usage to store result: >139.7 Mb
Dividing work into 128 chunks...
Total:9 s
vp <- fitExtractVarPartModel(vm_vp, f, d_all$ss, BPPARAM = bpparam())Dividing work into 128 chunks...
Total:9 s
saveRDS(vp, "input/vp.rds")
plotVarPart(sortCols(vp))f_ne <- ~ (1 | APOE) + (1 | siRNA) + (1 | stimulus) + (1 | line)
vm_vp_ne <- voomWithDreamWeights(d_all$d, f_ne, d_all$ss, BPPARAM = bpparam())Memory usage to store result: >139.7 Mb
Dividing work into 128 chunks...
Total:8 s
vp_ne <- fitExtractVarPartModel(vm_vp, f_ne, d_all$ss, BPPARAM = bpparam())Dividing work into 128 chunks...
Total:13 s
saveRDS(vp_ne, "input/vp_ne.rds")
plotVarPart(sortCols(vp_ne))Run Dream
run_dream <- function(ss, f, contrasts, cm, genes, grp,
norm_fun = calcNormFactors, fname = NULL) {
ss <- mutate(ss, across(where(is.factor), fct_drop))
if (!is.null(fname) && file.exists(fname)) return(read_rds(fname))
dgelist <- make_dgelist(ss, cm, genes, grp, norm_fun)
d <- dgelist[["d"]]
# Should we be using cpm(dgelist[["d"]]) instead of dgelist[["d"]]?
# transform count outcome variable to a continuous normally distributed
# outcome variable to run it in GLM
message("Making voom matrix")
vm <- voomWithDreamWeights(d, f, ss, BPPARAM = bpparam())
# make DREAM contrasts
if (!is.null(contrasts)) {
message("Making dream contrasts")
cd <- makeContrastsDream(f, ss, contrasts = contrasts)
# do the actual fitting of the linear model, running LMR one for each gene
message("fitting model using dream")
fit_raw = dream(vm, f, ss, cd, BPPARAM = bpparam())
} else {
fit_raw = dream(vm, f, ss, BPPARAM = bpparam())
cd = NULL
}
# eBayes - for each gene looks for genes with similar expression level and
# looks at its variance, looks for genes that varies much more than others
# and moderates their effect, red flag outliers would be color-coded, plot
# showing expression vs variance
message("applying bayesian shrinkage using Empirical Bayes")
fit = eBayes(fit_raw)
if (is.null(contrasts)) {
contrasts <- colnames(fit[["coefficients"]])
names(contrasts) <- contrasts
contrasts <- contrasts[contrasts != "rin"]
}
out <- list(inputs = list(ss = ss, cm = cm, f = f, contrasts = contrasts),
out = list(d = d, vm = vm, cd = cd, fit = fit, fit_raw = fit_raw,
matrix_dims = dgelist$matrix_dims))
if (!is.null(fname)) write_rds(out, fname)
return(out)
}
dr_genotype <-
run_dream(ss, ~ 0 + APOE + siRNA + stimulus + RIN + (1 | line_expt),
c("APOE" = "APOE44 - APOE33"),
cm, gencode45_genes, "APOE_siRNA_stimulus",
fname = "output/dr_APOE.rds")Making voom matrix
Memory usage to store result: >139.7 Mb
Dividing work into 128 chunks...
Total:6 s
Making dream contrasts
fitting model using dream
Dividing work into 128 chunks...
Total:374 s
applying bayesian shrinkage using Empirical Bayes
dr_siRNA <-
run_dream(ss, ~ 0 + siRNA + APOE + stimulus + RIN + (1 | line_expt),
c("siRNA_EIF2B2" = "siRNAEIF2B2- siRNASCR",
"siRNA_EIF2B3" = "siRNAEIF2B3- siRNASCR"),
cm, gencode45_genes, "APOE_siRNA_stimulus",
fname = "output/dr_siRNA.rds")Making voom matrix
Memory usage to store result: >139.7 Mb
Dividing work into 128 chunks...
Total:6 s
Making dream contrasts
fitting model using dream
Dividing work into 128 chunks...
Total:394 s
applying bayesian shrinkage using Empirical Bayes
dr_stimulus <-
run_dream(ss, ~ 0 + stimulus + siRNA + APOE + RIN + (1 | line_expt),
c("stimulus" = "stimulusTHPG- stimulusnone"),
cm, gencode45_genes, "APOE_siRNA_stimulus",
fname = "output/dr_stim.rds")Making voom matrix
Memory usage to store result: >139.7 Mb
Dividing work into 128 chunks...
Total:6 s
Making dream contrasts
fitting model using dream
Dividing work into 128 chunks...
Total:387 s
applying bayesian shrinkage using Empirical Bayes
dr_APOE_siRNA <-
run_dream(ss, ~ 0 + APOE_siRNA + stimulus + RIN + (1 | line_expt),
c("APOE_EIF2B2" = "APOE_siRNA44:EIF2B2 - APOE_siRNA33:EIF2B2",
"APOE_EIF2B3" = "APOE_siRNA44:EIF2B3 - APOE_siRNA33:EIF2B3"),
cm, gencode45_genes, "APOE_siRNA_stimulus",
fname = "output/dr_APOExsiRNA.rds")Making voom matrix
Memory usage to store result: >139.7 Mb
Dividing work into 128 chunks...
Total:6 s
Making dream contrasts
fitting model using dream
Dividing work into 128 chunks...
Total:395 s
applying bayesian shrinkage using Empirical Bayes
dr_APOE_siRNA_stim <-
run_dream(ss, ~ 0 + APOE_siRNA_stimulus + RIN + (1 | line_expt),
c("APOE_EIF2B2_THPG" = "APOE_siRNA_stimulus44:EIF2B2:THPG - APOE_siRNA_stimulus33:EIF2B2:THPG",
"APOE_EIF2B3_THPG" = "APOE_siRNA_stimulus44:EIF2B3:THPG - APOE_siRNA_stimulus33:EIF2B3:THPG"),
cm, gencode45_genes, "APOE_siRNA_stimulus",
fname = "output/dr_APOExsiRNAxstim.rds")Making voom matrix
Memory usage to store result: >139.7 Mb
Dividing work into 128 chunks...
Total:7 s
Making dream contrasts
fitting model using dream
Dividing work into 128 chunks...
Total:402 s
applying bayesian shrinkage using Empirical Bayes
plotSA(dr_genotype$out$fit)plotSA(dr_siRNA$out$fit)plotSA(dr_stimulus$out$fit)plotSA(dr_APOE_siRNA$out$fit)plotSA(dr_APOE_siRNA_stim$out$fit)GSEA
dgea <- function(fit, coef, genenames, padj.cutoff = 1,
out.dir = "output", write.xlsx = TRUE) {
if (is_tibble(fit)) {
res_unfilt <- fit |>
filter(contrast == coef)
} else {
res_initial <- topTable(fit, coef, n = Inf)
if (nrow(res_initial) == 0) return(NULL)
res_unfilt <- res_initial |>
rownames_to_column("gene_id") |>
left_join(genenames, by = "gene_id")
}
res <- res_unfilt |>
select(gene_id, gene_name, everything()) |>
filter(adj.P.Val <= padj.cutoff) |>
arrange(P.Value)
if (write.xlsx) res |> write_xlsx(str_glue("{out.dir}/{coef}.dgea.xlsx"))
return(res)
}
# Run GSEA on all of the contrasts within a DGE analysis
gsea <- function(fit, coef, genenames, pathways, padj.cutoff = 1,
in.dir = "input", out.dir = "output", write.xlsx = TRUE,
perms = 1000) {
stats <- dgea(fit, coef, genenames)
if (is.null(stats)) return(NULL)
stats <- stats |>
arrange(P.Value) |>
distinct(gene_name, .keep_all = TRUE) |>
arrange(desc(t))
stats <- setNames(stats$t, stats$gene_name)
res <- gmtPathways(str_glue("{in.dir}/{pathways}.symbols.gmt")) |>
fgsea(stats, eps = 0, BPPARAM = bpparam(), nPermSimple = perms) |>
as_tibble() |> # because fgsea returns a data.table object!
filter(padj <= padj.cutoff) |>
arrange(pval)
if (write.xlsx) res |> write_xlsx(str_glue("{out.dir}/{coef}.{pathways}.gsea.xlsx"))
return(res)
}
gsea_all <- function(results, genesets, gn = genenames,
padj.cutoff = 1, in.dir = "input", out.dir = "output",
fname = NULL, write.xlsx = TRUE, perms = 1000) {
if (!is.null(fname) & file.exists(file.path(out.dir, fname))) {
return(read_rds(file.path(out.dir, fname)))
} else if (!is.null(fname) & file.exists(fname)) {
return(read_rds(fname))
}
gsea_allpaths <- function(geneset, contrast, fit, genenames, padj.cutoff,
in.dir, out.dir, write.xls, perms) {
gsea(fit, contrast, gn, geneset, padj.cutoff, in.dir, out.dir,
write.xlsx, perms) |>
arrange(pval) |>
mutate(geneset = geneset, contrast = contrast,
pathway = if_else(str_detect(geneset, "^h\\.all\\.v"),
str_remove(str_replace_all(pathway, "_", " "),
"HALLMARK"),
pathway))
}
gsea_res <- expand_grid(x = genesets, y = names(results$inputs$contrasts)) |>
future_pmap_dfr(\(x, y) gsea_allpaths(x, y, results$out$fit, gn, padj.cutoff,
in.dir, out.dir, write.xlsx, perms),
.options = furrr_options(seed = TRUE))
out <- results
out$out$gsea <- gsea_res
if (!is.null(fname) & !str_detect(fname, "/")) {
write_rds(out, file.path(out.dir, fname))
} else if (!is.null(fname)) {
write_rds(out, fname)
}
out
}
# Specify which geneset files
gset = c("c2andc5.all.v2024.1.Hs", "h.all.v2024.1.Hs")
dr_genotype_gsea <- gsea_all(dr_genotype, gset, padj.cutoff = 1,
fname = "gsea_genotype.rds",
out.dir = "output",
perms = 10000)
dr_siRNA_gsea <- gsea_all(dr_siRNA, gset, padj.cutoff = 1,
fname = "gsea_siRNA.rds",
out.dir = "output",
perms = 10000)
dr_stimulus_gsea <- gsea_all(dr_stimulus, gset, padj.cutoff = 1,
fname = "gsea_stimulus.rds",
out.dir = "output",
perms = 10000)
dr_APOE_siRNA_gsea <- gsea_all(dr_APOE_siRNA, gset, padj.cutoff = 1,
fname = "gsea_APOE_siRNA.rds",
out.dir = "output",
perms = 10000)
dr_APOE_siRNA_stim_gsea <- gsea_all(dr_APOE_siRNA_stim, gset, padj.cutoff = 1,
fname = "gsea_APOE_siRNA_stim.rds",
out.dir = "output",
perms = 10000)
# alternative if you want the shiny without running GSEA
make_tables <- \(result, out.dir) {
fit_in <- result[["out"]][["fit"]]
contrasts <- names(result[["inputs"]][["contrasts"]])
fit_out <- map(contrasts, \(ctrst) {
dgea(fit_in, ctrst, genenames, out.dir = out.dir) |>
arrange(P.Value) |>
dgea_label()
})
names(fit_out) <- contrasts
fit_out
}Differential gene expression and gene set enrichment analyses
#summary table of degs across contrasts
summary(decideTests(dr_genotype$out$fit)) APOE APOE33 APOE44 siRNAEIF2B2 siRNAEIF2B3 stimulusTHPG RIN
Down 79 3063 3212 28 3 5740 2954
NotSig 14755 1675 1536 14851 14888 3369 9326
Up 57 10153 10143 12 0 5782 2611
#summary table of degs across contrasts
summary(decideTests(dr_siRNA$out$fit)) siRNA_EIF2B2 siRNA_EIF2B3 siRNASCR siRNAEIF2B2 siRNAEIF2B3 APOE44
Down 28 3 3063 3037 3023 79
NotSig 14851 14888 1675 1727 1704 14755
Up 12 0 10153 10127 10164 57
stimulusTHPG RIN
Down 5740 2954
NotSig 3369 9326
Up 5782 2611
#summary table of degs across contrasts
summary(decideTests(dr_stimulus$out$fit)) stimulus stimulusnone stimulusTHPG siRNAEIF2B2 siRNAEIF2B3 APOE44 RIN
Down 5740 3063 2865 28 3 79 2954
NotSig 3369 1675 2018 14851 14888 14755 9326
Up 5782 10153 10008 12 0 57 2611
#summary table of degs across contrasts
summary(decideTests(dr_APOE_siRNA$out$fit)) APOE_EIF2B2 APOE_EIF2B3 APOE_siRNA33:SCR APOE_siRNA33:EIF2B2
Down 52 1 3037 2964
NotSig 14813 14889 1716 1806
Up 26 1 10138 10121
APOE_siRNA33:EIF2B3 APOE_siRNA44:SCR APOE_siRNA44:EIF2B2
Down 2967 3133 3150
NotSig 1781 1635 1659
Up 10143 10123 10082
APOE_siRNA44:EIF2B3 stimulusTHPG RIN
Down 3133 5736 2945
NotSig 1606 3381 9337
Up 10152 5774 2609
#summary table of degs across contrasts
summary(decideTests(dr_APOE_siRNA_stim$out$fit)) APOE_EIF2B2_THPG APOE_EIF2B3_THPG APOE_siRNA_stimulus33:SCR:none
Down 551 384 2926
NotSig 14029 14303 1848
Up 311 204 10117
APOE_siRNA_stimulus33:SCR:THPG APOE_siRNA_stimulus33:EIF2B2:none
Down 2658 2820
NotSig 2264 1945
Up 9969 10126
APOE_siRNA_stimulus33:EIF2B2:THPG APOE_siRNA_stimulus33:EIF2B3:none
Down 2598 2753
NotSig 2359 2069
Up 9934 10069
APOE_siRNA_stimulus33:EIF2B3:THPG APOE_siRNA_stimulus44:SCR:none
Down 2694 3008
NotSig 2203 1809
Up 9994 10074
APOE_siRNA_stimulus44:SCR:THPG APOE_siRNA_stimulus44:EIF2B2:none
Down 2789 3027
NotSig 2218 1799
Up 9884 10065
APOE_siRNA_stimulus44:EIF2B2:THPG APOE_siRNA_stimulus44:EIF2B3:none
Down 2869 2951
NotSig 2194 1843
Up 9828 10097
APOE_siRNA_stimulus44:EIF2B3:THPG RIN
Down 2958 2961
NotSig 2055 9348
Up 9878 2582
res <- list("APOE" = dr_genotype_gsea,
"siRNA" = dr_siRNA_gsea,
"stimulus" = dr_stimulus_gsea,
"APOE_siRNA" = dr_APOE_siRNA_gsea,
"APOE_siRNA_stimulus" = dr_APOE_siRNA_stim_gsea) |>
map(\(result) {
fit_in <- result[["out"]][["fit"]]
contrasts <- names(result[["inputs"]][["contrasts"]])
fit_out <- map(contrasts, \(ctrst) {
dgea(fit_in, ctrst, genenames, write.xlsx = FALSE) |>
arrange(P.Value) |>
dgea_label()
})
names(fit_out) <- contrasts
result[["out"]][["fit"]] <- fit_out
if ("cd" %in% names(result[["out"]]) &&
"gsea" %in% names(result[["out"]])) {
result[["out"]] <- result[["out"]][c("fit", "gsea", "cd")]
} else if ("gsea" %in% names(result[["out"]])) {
result[["out"]] <- result[["out"]][c("fit", "gsea")]
} else if ("cd" %in% names(result[["out"]])) {
result[["out"]] <- result[["out"]][c("fit", "cd")]
} else {
result[["out"]] <- result[["out"]][c("fit")]
}
return(result)
})
shiny_title <- "B2 B3 KD APOE pop lines RNASeq"
write_rds(shiny_title, "shiny/title.rds")
save(res, genenames, mds, file = "results.rda")
R.utils::createLink("shiny/results.rda", "results.rda")[1] "shiny/results.rda"
Print environment
sessioninfo::session_info()─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.2.3 (2023-03-15)
os Ubuntu 22.04.4 LTS
system x86_64, linux-gnu
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz Etc/UTC
date 2025-01-14
pandoc 3.2.1 @ /sc/arion/work/goateomics/conda/envs/dreamcode/bin/ (via rmarkdown)
─ Packages ───────────────────────────────────────────────────────────────────
package * version date (UTC) lib source
abind 1.4-5 2016-07-21 [1] CRAN (R 4.2.3)
aod 1.3.3 2023-12-13 [1] CRAN (R 4.2.3)
backports 1.5.0 2024-05-23 [1] CRAN (R 4.2.3)
Biobase 2.58.0 2022-11-01 [1] Bioconductor
BiocGenerics 0.44.0 2022-11-01 [1] Bioconductor
BiocIO 1.8.0 2022-11-01 [1] Bioconductor
BiocParallel * 1.32.5 2022-12-23 [1] Bioconductor
Biostrings 2.66.0 2022-11-01 [1] Bioconductor
bit 4.0.5 2022-11-15 [1] CRAN (R 4.2.3)
bit64 4.0.5 2020-08-30 [1] CRAN (R 4.2.3)
bitops 1.0-7 2021-04-24 [1] CRAN (R 4.2.3)
boot 1.3-30 2024-02-26 [1] CRAN (R 4.2.3)
broom 1.0.6 2024-05-17 [1] CRAN (R 4.2.3)
car 3.1-2 2023-03-30 [1] CRAN (R 4.2.3)
carData 3.0-5 2022-01-06 [1] CRAN (R 4.2.3)
caTools 1.18.2 2021-03-28 [1] CRAN (R 4.2.3)
cellranger 1.1.0 2016-07-27 [1] CRAN (R 4.2.3)
cli 3.6.3 2024-06-21 [1] CRAN (R 4.2.3)
codetools 0.2-20 2024-03-31 [1] CRAN (R 4.2.3)
colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.2.3)
cowplot 1.1.3 2024-01-22 [1] CRAN (R 4.2.3)
crayon 1.5.3 2024-06-20 [1] CRAN (R 4.2.3)
crosstalk 1.2.1 2023-11-23 [1] CRAN (R 4.2.3)
data.table 1.15.4 2024-03-30 [1] CRAN (R 4.2.3)
DelayedArray 0.24.0 2022-11-01 [1] Bioconductor
digest 0.6.36 2024-06-23 [1] CRAN (R 4.2.3)
doParallel 1.0.17 2022-02-07 [1] CRAN (R 4.2.3)
dplyr * 1.1.4 2023-11-17 [1] CRAN (R 4.2.3)
edgeR * 3.40.0 2022-11-01 [1] Bioconductor
evaluate 0.24.0 2024-06-10 [1] CRAN (R 4.2.3)
fansi 1.0.6 2023-12-08 [1] CRAN (R 4.2.3)
farver 2.1.2 2024-05-13 [1] CRAN (R 4.2.3)
fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.2.3)
fastmatch 1.1-4 2023-08-18 [1] CRAN (R 4.2.3)
fgsea * 1.24.0 2022-11-01 [1] Bioconductor
forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.2.3)
foreach 1.5.2 2022-02-02 [1] CRAN (R 4.2.3)
furrr * 0.3.1 2022-08-15 [1] CRAN (R 4.2.3)
future * 1.33.2 2024-03-26 [1] CRAN (R 4.2.3)
generics 0.1.3 2022-07-05 [1] CRAN (R 4.2.3)
GenomeInfoDb 1.34.9 2023-02-02 [1] Bioconductor
GenomeInfoDbData 1.2.9 2024-08-29 [1] Bioconductor
GenomicAlignments 1.34.0 2022-11-01 [1] Bioconductor
GenomicRanges 1.50.0 2022-11-01 [1] Bioconductor
ggplot2 * 3.5.1 2024-04-23 [1] CRAN (R 4.2.3)
ggpubr * 0.6.0 2023-02-10 [1] CRAN (R 4.2.3)
ggrepel * 0.9.5 2024-01-10 [1] CRAN (R 4.2.3)
ggsignif 0.6.4 2022-10-13 [1] CRAN (R 4.2.3)
globals 0.16.3 2024-03-08 [1] CRAN (R 4.2.3)
glue 1.7.0 2024-01-09 [1] CRAN (R 4.2.3)
gplots 3.1.3.1 2024-02-02 [1] CRAN (R 4.2.3)
gtable 0.3.5 2024-04-22 [1] CRAN (R 4.2.3)
gtools 3.9.5 2023-11-20 [1] CRAN (R 4.2.3)
hms 1.1.3 2023-03-21 [1] CRAN (R 4.2.3)
htmltools 0.5.8.1 2024-04-04 [1] CRAN (R 4.2.3)
htmlwidgets 1.6.4 2023-12-06 [1] CRAN (R 4.2.3)
httr 1.4.7 2023-08-15 [1] CRAN (R 4.2.3)
IRanges 2.32.0 2022-11-01 [1] Bioconductor
iterators 1.0.14 2022-02-05 [1] CRAN (R 4.2.3)
jsonlite 1.8.8 2023-12-04 [1] CRAN (R 4.2.3)
KernSmooth 2.23-24 2024-05-17 [1] CRAN (R 4.2.3)
knitr 1.47 2024-05-29 [1] CRAN (R 4.2.3)
labeling 0.4.3 2023-08-29 [1] CRAN (R 4.2.3)
lattice 0.22-6 2024-03-20 [1] CRAN (R 4.2.3)
lazyeval 0.2.2 2019-03-15 [1] CRAN (R 4.2.3)
lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.2.3)
limma * 3.54.0 2022-11-01 [1] Bioconductor
listenv 0.9.1 2024-01-29 [1] CRAN (R 4.2.3)
lme4 1.1-35.3 2024-04-16 [1] CRAN (R 4.2.3)
lmerTest 3.1-3 2020-10-23 [1] CRAN (R 4.2.3)
locfit 1.5-9.9 2024-03-01 [1] CRAN (R 4.2.3)
magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.2.3)
MASS 7.3-60.0.1 2024-01-13 [1] CRAN (R 4.2.3)
Matrix 1.6-5 2024-01-11 [1] CRAN (R 4.2.3)
MatrixGenerics 1.10.0 2022-11-01 [1] Bioconductor
matrixStats 1.3.0 2024-04-11 [1] CRAN (R 4.2.3)
minqa 1.2.7 2024-05-20 [1] CRAN (R 4.2.3)
munsell 0.5.1 2024-04-01 [1] CRAN (R 4.2.3)
nlme 3.1-165 2024-06-06 [1] CRAN (R 4.2.3)
nloptr 2.0.3 2022-05-26 [1] CRAN (R 4.2.3)
numDeriv 2016.8-1.1 2019-06-06 [1] CRAN (R 4.2.3)
parallelly 1.37.1 2024-02-29 [1] CRAN (R 4.2.3)
pbkrtest 0.5.2 2023-01-19 [1] CRAN (R 4.2.3)
pillar 1.9.0 2023-03-22 [1] CRAN (R 4.2.3)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.2.3)
plotly * 4.10.4 2024-01-13 [1] CRAN (R 4.2.3)
plyr 1.8.9 2023-10-02 [1] CRAN (R 4.2.3)
prettyunits 1.2.0 2023-09-24 [1] CRAN (R 4.2.3)
progress 1.2.3 2023-12-06 [1] CRAN (R 4.2.3)
purrr * 1.0.2 2023-08-10 [1] CRAN (R 4.2.3)
R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.2.3)
R.oo 1.26.0 2024-01-24 [1] CRAN (R 4.2.3)
R.utils 2.12.3 2023-11-18 [1] CRAN (R 4.2.3)
R6 2.5.1 2021-08-19 [1] CRAN (R 4.2.3)
rbibutils 2.2.16 2023-10-25 [1] CRAN (R 4.2.3)
Rcpp 1.0.12 2024-01-09 [1] CRAN (R 4.2.3)
RCurl 1.98-1.14 2024-01-09 [1] CRAN (R 4.2.3)
Rdpack 2.6 2023-11-08 [1] CRAN (R 4.2.3)
readr * 2.1.5 2024-01-10 [1] CRAN (R 4.2.3)
readxl * 1.4.3 2023-07-06 [1] CRAN (R 4.2.3)
reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.2.3)
restfulr 0.0.15 2022-06-16 [1] CRAN (R 4.2.3)
RhpcBLASctl 0.23-42 2023-02-11 [1] CRAN (R 4.2.3)
rjson 0.2.21 2022-01-09 [1] CRAN (R 4.2.3)
rlang 1.1.4 2024-06-04 [1] CRAN (R 4.2.3)
rmarkdown 2.27 2024-05-17 [1] CRAN (R 4.2.3)
Rsamtools 2.14.0 2022-11-01 [1] Bioconductor
rstatix 0.7.2 2023-02-01 [1] CRAN (R 4.2.3)
rstudioapi 0.16.0 2024-03-24 [1] CRAN (R 4.2.3)
rtracklayer 1.58.0 2022-11-01 [1] Bioconductor
S4Vectors 0.36.0 2022-11-01 [1] Bioconductor
scales 1.3.0 2023-11-28 [1] CRAN (R 4.2.3)
sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.2.3)
stringi 1.8.4 2024-05-06 [1] CRAN (R 4.2.3)
stringr * 1.5.1 2023-11-14 [1] CRAN (R 4.2.3)
SummarizedExperiment 1.28.0 2022-11-01 [1] Bioconductor
tibble * 3.2.1 2023-03-20 [1] CRAN (R 4.2.3)
tidyr * 1.3.1 2024-01-24 [1] CRAN (R 4.2.3)
tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.2.3)
tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.2.3)
utf8 1.2.4 2023-10-22 [1] CRAN (R 4.2.3)
variancePartition * 1.28.0 2022-11-01 [1] Bioconductor
vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.2.3)
viridisLite 0.4.2 2023-05-02 [1] CRAN (R 4.2.3)
vroom 1.6.5 2023-12-05 [1] CRAN (R 4.2.3)
withr 3.0.0 2024-01-16 [1] CRAN (R 4.2.3)
writexl * 1.5.0 2024-02-09 [1] CRAN (R 4.2.3)
xfun 0.45 2024-06-16 [1] CRAN (R 4.2.3)
XML 3.99-0.17 2024-06-25 [1] CRAN (R 4.2.3)
XVector 0.38.0 2022-11-01 [1] Bioconductor
yaml 2.3.8 2023-12-11 [1] CRAN (R 4.2.3)
zlibbioc 1.44.0 2022-11-01 [1] Bioconductor
[1] /sc/arion/work/goateomics/conda/envs/dreamcode/lib/R/library
──────────────────────────────────────────────────────────────────────────────